Finite-size effects at first-order isotropic-to-nematic transitions 
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We present simulation data of first-order isotropic-to-nematic transitions in lattice models of 
liquid crystals and locate the thermodynamic limit inverse transition temperature e^o via finite-size 
scaling. We observe that the inverse temperature of the specific heat maximum can be consistently 
extrapolated to toe assuming the usual a/L'^ dependence, with L the system size, d the lattice 
dimension and proportionality constant a. We also investigate the quantity eL,fe, the finite-size 
inverse temperature where k is the ratio of weights of the isotropic to nematic phase. For an 
optimal value k — fcopt, ^L,k versus L converges to Coo much faster than a/L'^ , providing an economic 
alternative to locate the transition. Moreover, we find that a ~ Inkopt/ Coa, with Coo the latent 
heat density. This suggests that liquid crystals at first-order IN transitions scale approximately as 
g-state Potts models with q ~ fcopt. 



PACS numbers: 05.70.Fh, 75.10.Hk, 64.60.-i, 64.70.mf 

I. INTRODUCTION 

The investigation of the isotropic-to-nematic (IN) tran- 
sition in liquid crystals via computer simulation is long 
established. Decades ago Lebwohl and Lasher (LL) in- 
troduced a simple lattice model, the LL model, to study 
this transition At each site z of a cubic lattice they 
attached a three-dimensional unit vector di (spin) inter- 
acting with its nearest-neighbors via 



(1) 



<i]> 



where p — 2 and with a factor l/fc^T absorbed into the 
coupling constant e > 0, with fcs the Boltzmann constant 
and T the temperature. Despite its simplicity, the LL 
model captures certain aspects of liquid crystal phase be- 
haviour remarkably welL and has consequently received 
considerable attention pi, l3| . 

A common problem in locating the IN transition via 
simulation is the issue of finite system size. Phase tran- 
sitions are defined in the thermodynamic limit, whereas 
simulations always deal with finite particle numbers. In 
order to estimate the thermodynamic limit transition 
point, it is typical to perform a number of simulations 
for different system sizes and to subsequently extrapo- 
late the results following some finite-size scaling (FSS) 
procedure. Which procedure to use depends on the type 
of transition, i.e. whether it is continuous or first-order. 
In three spatial dimensions, the IN transition is typically 
first-order; in two dimensions both continuous jj, [1, [6| 
and first-order 0,[1] IN transitions can occur, depending 
on the details of the interactions [1, [l^, [ll| . In this paper 
we focus on the first-order case. 

The literature on FSS at first-order transitions is quite 
extensive, for a review see [l^ . the majority of which 
deals exclusively with the Potts model [ij]. An impor- 
tant result is that the "apparent" transition inverse tem- 
perature ei.cv, obtained in a finite system of size L, is 



shifted from the thermodynamic limit value too as [T3, ITsj 



f-L,CV 
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with proportionality constant a > 0. Here eL,cv is the 
inverse temperature where the specific heat in a finite 
system of size L attains its maximum, d is the spatial di- 
mension of the lattice and L denotes the linear extension 
of the simulation box, generally square or cubic, with 
periodic boundary conditions. 

We emphasize that Eq.([21) was derived for the Potts 
model where the proportionality constant is known to be 



apotts = ln(7/£o 



(3) 



with q the number of Potts states and C^n the latent 
heat density in the thermodynamic limit [T4 . Il5| . Inter- 
estingly, simulations of the LL model have shown that 
the functional form of Eq. ^ also works well for IN tran- 
sitions [l^ ■ That is, meaningful extrapolations of 
^L,cv can be performed, although the significance of a is 
not obvious. It certainly cannot be related to the number 
of spin states, i.e. conform to Eq.Q, since the LL model 
is a continuous spin model, in contrast to the discrete 
spin variables of the Potts model ■ 

In any case, based on the success of Eq.([2]) in describing 
finite-size effects in the LL model, it could be hoped that 
other scaling relations, originally derived for the Potts 
model, also remain valid. Of particular interest is the 
result of Borgs and Kotecky, who showed that for the 
Potts model exponentially decaying finite-size effects are 
also possible [isl fiol . [20| . The obvious advantage of expo- 
nential decay is that too is approached much faster with 
increasing L, compared to the power law decay of Eq.(I2). 
This means that moderate system sizes may suffice to 
locate the transition, thereby saving valuable computer 
time. As liquid crystal phase transitions are in any case 
expensive to simulate, such a gain in efficiency would cer- 
tainly be highly desirable. 

We will show in this paper that it is indeed possible to 
locate first-order IN transitions from finite-size Simula- 
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tion data with shifts that vanish much faster than \/L'^. 
This is possible by considering e^^fc, the inverse tempera- 
ture at which the "ratio-of- weights" of the isotropic and 
the nematic phases is equal to a value k. This ratio- 
of-wcights is obtained from the order parameter distri- 
bution PL.e{X), defined as the probability to observe an 
order parameter A, when simulating a system of size L 
at inverse temperature e. In the vicinity of the IN tran- 
sition the distribution becomes bimodal, with one peak 
corresponding to the isotropic phase and the other to the 
nematic phase. The ratio-of-weights is simply the ratio of 
the peak areas. Provided k is chosen optimally eL,k ap- 
proaches eoo extremely rapidly as L increases, yielding an 
economic alternative over Eq.([2l). A prerequisite is that 
the transition must be strong enough first-order for the 
ratio-of-weights to be meaningfully calculated [21(|. For 
this reason, we do not consider the original LL model, as 
the transition is extremely weak here, but a variation of 
it. 

In this paper we firstly provide the details of the mod- 
ified LL model in Section |TT1 together with a description 
of the simulation method that was used to obtain the 
order parameter distribution. Next, we measure eoo us- 
ing the "standard" approach of extrapolating €l^cv via 
Eq.(I21), as well as using the "new" approach based on 
^L,k- In particular we demonstrate how to locate the op- 
timal value fcopt, along which finite-size effects are min- 
imal. As expected, both approaches are in good agree- 
ment, with the essential difference that e^^k converges to 
Eoo already for very small systems. This fast convergence 
property was observed at all transitions studied by us, ir- 
respective of space and spin dimension. We also consider 
the finite-size scaling of the latent heat density and show 
that, for IN transitions, fcopt becomes the "analogue" of 
the number of Potts states q. Finally we present a sum- 
mary of our findings in Section IIVI 



haviour is observed in the LL model: even in simulation 
boxes of L = 70 lattice spacings the peaks strongly over- 
lap and the logarithm of the peak height, measured with 
resp ect to the minimum in between, is less than 2 /cbT 
[l7j . Since the peaks overlap one never truly sees pure 
phases, which complicates the analysis. In order to yield 
reasonable results we require in this paper that the peaks 
in PL,e(A) be well-separated. More precisely, it must be 
possible to assign a "cut-off" separating the peaks, on 
which the final results may not sensitively depend. For 
this reason we do not consider the original LL model but 
rather a generalization of it, where the exponent p of 
Eq.(IT|) exceeds the LL value. We expect this will lead to 
a much stronger first-order IN transition [l^, [ll[ , so 
distributions will display non-overlapping peaks already 
in moderately sized systems. In fact by using a large ex- 
ponent p in Eq.((T]) strong first-order IN transitions rnay 
be realized even in purely two-dimensional systems 0, H] . 
Hence, the model that we consider is just the LL model 
of Eq.([T|) but with p > 2. Note the absolute value | • | 
such that the system is invariant under inversion of the 
spin orientation. We thus impose the symmetry of liquid 
crystals although we believe that our results also apply 
to magnetic systems. 

Note that the use of a large exponent p in Eq. ([1]) may 
also yield a better description of experiments on con- 
fined liquid crystals. The latter systems are quasi two- 
dimensional. If one studies the LL model in two dimen- 
sions, i.e. with p — 2, and three-dimensional spins, a true 
phase transition appears to be absent (23i |. In contrast, 
experiments clearly reveal that transitions do occur. In 
fact, these transitions appear to be of the IN type and 
are quite strong, as manifested by pronounced coexis- 
tence between isotropic and nematic domains [24]. Such 
behaviour cannot be reproduced easily with the standard 
LL model, but it can be using the modified version con- 
sidered in this work, with a sufficiently large exponent p. 



II. MODEL AND SIMULATION METHOD 
A. Modified LL model 

In order to study finite-size effects at phase transitions, 
simulation data of high statistical quality are essential. 
This sets a limit on the complexity of the models that 
can be handled as well as on the system size. For our 
purposes already the simple LL model is too demand- 
ing, the problem being that the IN transition in this 
model is extremely weak. Generally, in computer sim- 
ulations, first-order phase transitions are identified by 
measuriiig the probability distribution of the order pa- 
rameter [22| . At the transition, this distribution displays 
two peaks: one corresponding to the isotropic phase and 
the other to the nematic phase. In the thermodynamic 
limit, the peaks become sharper, and ultimately a distri- 
bution of two (5-functions is obtained. In finite systems, 
however, the peaks are broad and possibly overlapping, 
especially when the transition is very weak. Such be- 



B. Transition matrix Wang-Landau sampling 

Following earlier work on the LL model [l^, our 
simulations are based on the order-parameter distribu- 
tion. We use the energy E of Eq. ([T]) as order parameter 
and aim to measure PL,e{E) as accurately as possible. 
Recall that PL,g{E) is the probability to observe energy 
E, in a system of size L, at inverse temperature e. De- 
pending on the case of interest, the simulations are per- 
formed on square or cubic lattices of linear size L, using 
periodic boundary conditions. 

In order to obtain Pl^^{E) we use Wang-Landau (WL) 
sampling [2^ [2^ additionally optimized by recording 
some elements of the transition matrix (TM) 27, 2^. 
The aim of WL sampling is to perform a random walk 
in energy space, such that all energies are visited equally 
often. To this end, we use single spin dynamics, whereby 
one of the spins is chosen randomly and given a new ran- 
dom orientation. The new state is accepted with proba- 
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bility 



p{Ei Ej) = mill 



9{Ej 



1 



(4) 



with Ej and Ej the energies of the initial and final states 
respectively and g{E) the density of states. The density 
of states is unknown beforehand and g{E) is initially set 
so g{E) = 1. Upon visiting any particular energy the 
corresponding density of states is multiplied by a modifi- 
cation factor / > 1. We also keep track of the histogram 
H{E), counting the number of times each energy E is 
visited. Once H{E) contains sufhcient information over 
the range of energy of interest, the modification factor 
/ is reduced and the energy histogram H{E) is reset to 
zero. These steps are repeated until / has become close 
to unity, after which changes in the density of states be- 
come negligible. The sought order parameter distribution 
is then obtained from Pl,^{E) cx g{E) exp{~eE). 

The above procedure is the standard WL algorithm, 
which works extremely well in many cases 1]. However, 
it has been noted [H, [2^ that the WL algorithm in its 
standard form reaches a limiting accuracy, after which 
the statistical quality of the data no longer improves, no 
matter how much additional computer time is invested. 
Hence, these authors also propose to measure the TM 
elements T{Ei Ej). These are defined as the number 
of times that, being in a state with energy Ei, a state 
with energy Ej is proposed, irrespective of whether the 
new state is accepted. From the TM elements one can 
estimate 



n{Ei Ej) 



T{Ei 



El 



EKTiEi-^EK) 



(5) 



which is the probability that being in state with energy 
Ej, a move to a state with energy Ej is proposed. This 
is related to the density of states via 



g{Ei) ^ n{Ej ^ Ei) 

g{Ej) n{Ei ^ Ej) ■ 



(6) 



Hence, by recording TM elements the density of states 
can also be constructed, the great advantage being that 
rejected moves also give useful information. 

To combine WL sampling with the TM method we 
somewhat follow (2^ . At the start of the simulation the 
density of states g{E) is set to unity, while the energy 
histogram H{E) and the TM elements are set to zero. 
We perform one WL iteration, i.e. accepting moves con- 
form Eq.Q, using a high modification factor In/ = 1. 
At each move both H{E) and the TM elements are up- 
dated. We continue to simulate until all bins in H{E) 
contain at least n entries over the chosen energy range. 
We then use the TM elements to construct a new density 
of states, which serves as the starting density of states for 
the next WL iteration. For the next iteration H{E) is re- 
set to zero, the modification factor is reduced to {\nf )/l 
but the TM elements remain untouched. These steps are 
repeated until In / w 0(10"^°), after which we store the 



corresponding density of states gp{E). This marks the 
end of the "prepare" stage. 

Next we proceed with the "collect" stage. The TM ele- 
ments are set to zero, whereas H{E) is no longer needed. 
During collection we sample according to Eq.Q using 
gp{E) as estimate for the density of states. However, 
only the TM elements and not gp{E) are further up- 
dated. As collection proceeds the accuracy of the TM el- 
ements increases indefinitely, as does the accuracy of the 
density of states obtained from them. The reason to have 
a separate "collect" stage is because during "prepare" de- 
tailed balance is not strictly obeyed, due to the initially 
large modification factor / [2^. For this reason, gp{E) 
could be biased and we are reluctant to perform finite- 
size scaling with it. 

During the "prepare" stage small values n w 10 to- 
gether with large values I « 5 — 10 can be used. This 
significantly speeds up the simulation and similar obser- 
vations have been made in other works ^28<] . Histograms 
were collected by discretizing the energy in bins of reso- 
lution A = 1 kpT. In order to avoid "boundary effects" 
during WL sampling states are counted as in [soj . To 
reduce memory consumption, only the nearest-neighbor 
elements T{Ej Ej ± A) of the TM along with the 
normalization T{Ej Ek) are stored. Since sin- 
gle spin dynamics are used, these are the dominant en- 
tries. Constructing the density of states using Eq.® 
and recursion is then a straightforward matter. If all 
TM elements were to be used, constructing the density 
of states becomes more complex while not yielding sig- 
nificantly higher accuracy [23], so this is not attempted 
here. The required computer time depends sensitively on 
the size of the system. For small systems consisting of 
~ 1000 spins, the "prepare" stage can be completed in 
as short a time as 15 minutes. For larger systems con- 
taining 10, 000 spins or more this can take more than one 
week. In these cases it is necessary to collect the density 
of states over a number of separate energy intervals with 
a single processor assigned to each interval. Such a par- 
allelization is trivially implemented. The "collect" phase 
typically lasts as long as the "prepare" phase except for 
very large systems, where it is found that it takes a much 
longer time to obtain an equivalently accurate density of 
states. 



III. RESULTS AND ANALYSIS 

We have performed extensive simulations of Eq.((T]) 
varying both the space and spin dimension, as well as 
the exponent p. More precisely, the following scenarios 
are considered: 

1. three-dimensional lattices, three-dimensional spins 
with p = 5 — 45, 

2. two-dimensional lattices, three-dimensional spins 
with p = 20 — 50 and 
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FIG. 1: Logarithm of Ph,t(E) using p = 10 in Eq.([T} for sys- 
tem sizes L = 10, 12, 15 (from top to bottom), cubic lattices, 
and three-dimensional spins. In each of the distributions e 
was tuned so the peaks are of equal height. The barrier A_F, 
here marked for the L = 10 system, is defined as the height of 
the peaks measured with respect to the minimum in between. 
The peak-to-peak distance Ap corresponds to the latent heat 
density. Note that we have plotted the distributions as a 
function of the negative energy density: the left peak thus 
corresponds to the isotropic phase and the right peak to the 
nematic phase. 



3. two-dimensional lattices, 
with p = 150 - 1000. 



two-dimensional spins 



On three-dimensional lattices, it is well accepted that the 
IN transition is first-order. The fact that the IN transi- 
tion can also be first-order in two dimensions is perhaps 
less well known. In this case, first-order transitions only 
appear pro vided the exponent -p of Eq.(IT]) is sufficiently 
large [ol. lid. flit. Hence, in two dimensions, one generally 
needs p ^ 2 in order to observe a first-order transition, 
and it is important to verify that such a transition is 
indeed taking place. If one additionally lowers the spin 
dimension from 3^2, even greater exponents p are re- 
quired. For this reason, the chosen p ranges vary signifi- 
cantly between the three scenarios. 



A. Determining the order of the transition 

In order to verify the presence of a first-order transition 
we use the scaling method of Lee and Kosterlitz 3l|, [l^l . 
Recall that the order parameter distribution becomes bi- 
modal in the vicinity of the IN transition; see Fig. [1] for 
an example. The idea of Lee and Kosterlitz is to monitor 
the peak heights AF of the logarithm of the order param- 
eter distribution, measured with respect to the minimum 
"in-between" the peaks. At a first-order transition AF 
corresponds to the formation of interfaces between co- 
existing isotropic and nematic domains j33j . In d spatial 
dimensions it is therefore expected that AF oc L'^~^, pro- 



FIG. 2: Variation of AF versus L''"^ for d = 3 dimensional 
lattices, three-dimensional spins, and the various values of p 
as indicated. 
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FIG. 3: Similar to Fig.[2]but for two-dimensional lattices and 
three-dimensional spins. 
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FIG. 4: Similar to Fig.[2]but for two-dimensional lattices and 
two-dimensional spins. 
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FIG. 5: Estimation of e^o via extrapolation of eL,cv using 
Eq.(l2]). Shown is eL,cv versus l/L'* using (a) p = 10, cu- 
bic lattices and three-dimensional spins, (b) p = 20, square 
lattices and three-dimensional spins and (c) p = 150, square 
lattices and two-dimensional spins. The open symbols are 
simulation data and the lines are fits to Eq.(l2]). 



viding a straightforward recipe to identify the transition 
type: a linear increase of AF versus L''"^ indicates that a 
first-order transition is taking place, with the slope yield- 
ing the interfacial tension 33], whereas for a continuous 
transition AF becomes independent of L, or vanishes 
altogether if no transition takes place at all in the ther- 
modynamic limit. In Fig. [2]AF is plotted for the purely 
three-dimensional case; the linear increase is clearly vis- 
ible, confirming the presence of a first-order transition. 
For two-dimensional lattices the results have been col- 
lected in Fig. [3] and Fig. [D Once again the presence of 
a first-order transition is confirmed. Note that on two- 
dimensional lattices the slopes of the lines correspond to 
line tensions. 



B. Extrapolation of eL,cv 

Next, we measure the thermodynamic limit inverse 
temperature e^o by means of extrapolating eL.cv via 
Eq.([2)- Recall that eL,cv is the finite-size inverse tem- 



{E') {EY 



(7) 



attains its maximum. Shown in Fig.[5fa) is eL,cv versus 
for the purely three-dimensional case using p — 
10 - results for different p are qualitatively similar and 
therefore not explicitly shown. In agreement with earlier 
simulations of the "original" LL model [TgI . , the data 
are well described by Eq.Q and from the fit eoo can be 
meaningfully obtained. The resulting fit parameters are 
collected in Table HI Repeating the same analysis for two- 
dimensional lattices yields similar results; some typical 
plots are shown in Fig. E^b) and (c) with the resulting fit 
parameters collected in Tables II and III. 



C. Extrapolation of e^^k 

We now arrive at the main result of this paper, namely 
the estimation of too by monitoring eL,k- Recall that e^ k 
is defined as the finite-size inverse temperature where the 
equality 



Wn/Wi = k 



(8) 



is obeyed, with Wn and Wi the areas under the nematic 
and isotropic peaks of the order parameter distribution 
respectively. No matter what value of k is used, provided 
it is positive and finite, we expect that linii^oo eL,k = 
Eoo- The reason is that in the thermodynamic limit a 
bimodal order parameter distribution survives only at 
Eoo and not anywhere else [33 |. Hence, keeping the area 
ratio fixed at some value of k whilst increasing L, e^ k 
will definitely approach ■ The rate of the convergence, 
however, does depend on k. Assuming that the prediction 
of Borgs and Kotecky for the Potts model also holds at 
first-order IN transitions, it should be possible to locate 
an optimal value /copt at which the convergence to eoo 
is fastest and hopefully faster than 1/L''. Therefore, we 
propose to manually inspect the convergence of eL,fc using 
several values of k. 

A prerequisite for numerically solving Eq.((8|) is that 
the transition must be sufficiently first-order in order for 
bimodal distributions PL,t{E) with well-separated peaks 
to appear. By this we mean that the barrier AF ^ defined 
in Fig. [U is large enough. The areas of the nematic and 
isotropic peaks may then be calculated using 

Wn^ PLAE)dE,Wi^ PL.eiE)dE, (9) 

where we remind the reader that the energy in our model 
is negative. The details of defining the "cut-off" energy 
Ec are somewhat arbitrary, but as states around Ec con- 
tribute exponentially little to the peak areas, the precise 
form does not matter 34] . In this work Ec is taken to be 
the average Ec = J EPl^E) dE, with Pl^E) obtained 



6 



TABLE L Properties of the IN transition of Eq.([T| for three-dimensional lattices and three-dimensional spins versus p. Listed 
are the fit parameters eoo,cv and a of Eq.([2]), the best estimate eoo,fc obtained from the convergence of eL.fc along fcopt, the 
logarithm of fcopt with uncertainty AA;, the latent heat density £00 and the ratio In fcopt /-Coo • 



p 




a 


^00 , A: 


In fcopt ± Afc 


Coo 


In fcopt / £00 


5 


1.3969 


6.62 


1.3970 ±0.001 


2.7 ±0.6 


0.358 ±0.010 


5.9 - 9.2 


8 


1.5207 


5.28 


1.5207 ±0.001 


5.0 ±0.2 


0.908 ± 0.005 


5.3 - 5.7 


10 


1.5862 


4.28 


1.5864 ±0.001 


5.3 ±0.5 


1.156 ±0.005 


4.2-5.0 


15 


1.7126 


3.94 


1.7126 ±0.001 


6.0 ±0.2 


1.523 ±0.002 


3.8-4.1 


20 


1.8063 


3.76 


1.8063 ±0.001 


6.4 ±0.3 


1.727 ±0.002 


3.5-3.9 


45 


2.0838 


3.55 


2.0838 ±0.001 


8.1 ± 1.0 


2.120 ±0.001 


3.3-4.3 



TABLE IL Similar to Table |T] but for two-dimensional lattices and three-dimensional spins. 



p 


£00, cv 


Q 




In fcopt ± Afc 




In /copt /Coo 


20 


2.7695 


5.18 


2.7698 ±0.001 


4.1 ±0.4 


0.7145 ±0.001 


5.2-6.3 


25 


2.8678 


4.84 


2.8679 ± 0.001 


4.5 ±0.2 


0.8900 ± 0.001 


4.8 - 5.3 


30 


2.9517 


4.69 


2.9517 ±0.001 


4.8 ±0.2 


1.0023 ±0.001 


4.6 - 5.0 


35 


3.0240 


4.58 


3.0241 ±0.001 


5.2 ±0.6 


1.0832 ±0.001 


4.2 - 5.4 


40 


3.0882 


4.52 


3.0882 ± 0.001 


5.2 ±0.2 


1.1432 ±0.001 


4.4 - 4.7 


45 


3.1456 


4.47 


3.1455 ±0.001 


5.2 ±0.5 


1.1936 ±0.001 


3.9-4.8 


50 


3.1976 


4.50 


3.1976 ±0.001 


5.6 ±0.5 


1.2320 ±0.001 


4.1 - 5.0 



at equal-height, i.e. as in Fig. [TJ Once has been set 
its value is kept fixed whilst solving Eq.®. 

For the purely three-dimensional case, the behaviour 
of eL,fc is shown in Fig. [SI Using a number of exponents 
p in Eq.((T]), we have plotted tL,k versus L for several val- 
ues of k. The data are consistent with the expectation 
that, regardless of fc, t^^k converges to a common value, 
corresponding to eoo- Note also that Coo is approached 
from above for large k, and from below for low k. Hence, 
we can indeed identify an optimal value fcopt along which 
finite-size effects are minimal. The optimum can be esti- 
mated by locating, for a pair of system sizes Li and Lj , 
the inverse temperature where for both system sizes 
the same ratio fcy of the peak areas is observed. By con- 
sidering all available pairs of system sizes, the average 
and root -mean-square fluctuation in e^j and k^j can be 
calculated, which then yield Coq and fcopt with uncertain- 
ties, shown in Table [J Although fcopt itself is not known 
very precisely, since Afc is quite large, very accurate esti- 
mates of £00 can still be obtained as this quantity is rather 
insensitive to the precise value of fcopt being used. This 
means that the series tL,k also provides a valid method 
for locating IN transitions. The corresponding estimates 
of 600 are in good agreement with those obtained via ex- 
trapolation of tL^cVi as inspection of the various tables 
indicates. The practical advantage of using t^^k with 
fc = fcopt is that the L-dependence is very weak, so much 
so that £00 is captured already in small systems. Similar 
findings were obtained using two-dimensional lattices, of 
which some typical plots are provided in Fig. [7] and Fig.[H] 
with the corresponding numerical estimates collected in 



Tables II and III. 

For non-optimal values fc ^ fcopt, we observe that the 
shift £00 — £L,fc oc l/i'', i.e. the shift vanishes as a power 
law in the inverse volume, similar to ^l.cv- At the opti- 
mal value fc = fcopt, finite-size effects in e^^^ are typically 
too small in order for a meaningful fit to be carried out. 
Hence, our data confirm Borgs and Kotecky in the sense 
that optimal estimators can be defined which converge 
onto £00 faster than 1/L'^\ whether the optimal conver- 
gence is indeed exponential requires more accurate data, 
which is currently beyond our reach [ssl ]. 

An alternative, but completely equivalent, method 
to investigate the convergence of e^^k is presented in 
[H, [36, 37, 38], albeit for the Potts model. The idea 
is to plot the area ratio Wn/Wj versus e for several sys- 
tem sizes. The resulting curves are expected to reveal 
an intersection point at the transition inverse tempera- 
ture; the value of the area ratio at the intersection then 
yields fcopt • For completeness we have prepared one such 
plot, see Fig. [SI The curves indeed intersect and give 
estimates of £00 and fcopt that are fully consistent with 
those reported in Table I. 



D. Latent heat density 

It appears that scaling relations derived for the Potts 
model also work remarkably well at IN transitions. In 
agreement with earlier simulations of the LL model [igI . 
iiZl, the validity of Eq.Q is confirmed additionally by 
us. Furthermore, our data suggest that an analogue of 
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TABLE III: Similar to Table |T] but for two-dimensional lattices and two-dimensional spins. 



V 


Eoo.CV 


a 


^OO , k 


In fcopt ± Afc 




In /<^opt/-^oo 


150 


2.063 


3.67 


2.0628 ± 0.0005 


1.9 ±0.8 


0.648 ± 0.005 


1.7 - 4.2 


1000 


2.486 


3.31 


2.4865 ± 0.0006 


4.8 ± 1.3 


1.270 + 0.005 


2.8 - 4.8 




FIG. 6: Variation of e^.fc versus L, for three-dimensional lat- 
tices and three-dimensional spins, using different exponents 
p as indicated. The symbols are simulation data, the lines 
serve to guide the eye. The central curves in each plot show 
eL,k using k = fcopt along which finite-size effects are minimal; 
also shown is eL,k using fc = fcopt — 5 Afc (lower curves), and 
fc = fcopt + 5 Afc (upper curves). The methods for locating fcopt 
and Afc are explained in the text and the resulting values, as 
well as the estimates of eao, are listed in Table U 



the Borgs and Kotecky prediction, namely that finite-size 
effects vanish faster than l/L"* at appropriate points, can 
be defined. In this case, it is needed to measure tL,k 
using the optimal value k = A:opt- In the Potts model 
it holds that fcopt = 9, where q is the number of Potts 
states. In other words, finite-size effects in the Potts 
model are minimized when the ratio of the peak areas in 
the order parameter distribution is held fixed at q. Based 
on our results, it seems reasonable to assume that scaling 
relations for the Potts model also hold at IN transitions, 
but with q replaced by fcopt- 

To test this assumption we consider the proportional- 
ity constant a from Eq.Q, which is given by Eq.(I3]) for 
the Potts model. If q can be replaced by fcopt , ol should 
correspond to In fcopt /Coo , where £oo is the latent heat 
density. The latter can be obtained independently from 
C'l, max = where CL,max is the maximum value 



FIG. 7: Similar to Fig. [6] but for two-dimensional lattices 
and three-dimensional spins. Numerical estimates are given 
in Table nil 




2.495 



2.490 



2.485 



2.480 



FIG. 8: Similar to Fig. (6] but for two-dimensional lattices 
and two-dimensional spins. Numerical estimates are given in 
Table Iml 
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of the specific heat in a finite system of size L [15|. Hence, 
we introduce the latent heat estimator 



1.5850 



1.5870 



FIG. 9; Plot of \nWN /Wi versus e for several system sizes 
L. This plot was generated using p = 10 in Eq.([l}, three- 
dimensional lattices and three-dimensional spins. 




100 



120 



FIG. 10: Finite-size variation of the latent heat density es- 
timators Api,,i and ApL,2. Results are shown for Eq.(IT| 
using (a) p = 10, three-dimensional lattices and three- 
dimensional spins, (b) p — 20, two-dimensional lattices and 
three-dimensional spins and (c) p = 150, two-dimensional lat- 
tices and two-dimensional spins. 



ApL,l = \/4CL^max/i^ 



(10) 



which should approach Coo as i ^ oo. Addition- 
ally, the latent heat density can be read-off directly, 
as the peak-to-peak distance in the energy distribution, 
marked Ap in Fig. [TJ Numerically this is expressed by 
Ml = '^{\E - {E)\)/L'^] plotting Ml versus e gives a 
maximum /S.pL,2, which in the limit L — > oo also ap- 
proaches Coo- Typical behaviour of A.pL,i is shown in 
Fig. [ini As expected, both latent heat estimators con- 
verge to a common value, which can be read-off reason- 
ably accurately; the resulting estimates of Coo are given 
in the various tables. Note also that Coo is approached 
from below in three dimensions, whereas in two dimen- 
sions, it is approached from above. If an appropriate 
number of two-dimensional lattice layers stacked on top 
of each other were simulated, it is likely that a cross-over 
regime could be found where ApL,i depends only weakly 
on L, as these systems are effectively in-between two and 
three dimensions. 

Having measured Coo, the ratio Infcopt/'Coo is eas- 
ily obtained, which may then be compared to a, 
see Tables I-III. The uncertainty is admittedly rather 
large, but within numerical precision, and the relation 
In fcopt /Coo ~ C( appears to hold. 



IV. SUMMARY 

In this paper we have presented simulation data of 
first-order isotropic-to-nematic transitions in lattice liq- 
uid crystals with continuous orientational degrees of free- 
dom for various space and spin dimensions. As with ear- 
lier simulations of this type [l^, [l3| , we find that the ex- 
trapolation of the finite-size inverse temperature of the 
specific heat maximum eL,cv can be consistently per- 
formed assuming a leading a/L'^ dependence, exactly as 
in the Potts model. Inspired by this result, we have inves- 
tigated an alternative approach to locate the transition 
inverse temperature using estimators eL,fe, defined as the 
finite-size inverse temperature where the ratio of peak ar- 
eas in the energy distribution is equal to k. In agreement 
with the Potts model, eL,k converges to Coo much faster 
than 1/L'^, provided an optimal value k = fcopt is used. 
Moreover, the ratio kopt/Coo, with Coo the latent heat 
density, is remarkably consistent with the proportional- 
ity constant a from the scaling of cl.cv- This leads us 
to conclude that finite-size scaling predictions originally 
proposed for first-order transitions in the Potts model re- 
main valid at first-order IN transitions too, but with the 
number of Potts states q replaced by fcopt ■ 

It is perhaps somewhat surprising that a continuous 
spin model at a first-order transition, such as the LL 
model, scales in the same way as the Potts model, which 
is, after all, a discrete spin model. In fact, Borgs and 
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Kotecky have remarked that the derivation of their scal- 
ing results cannot be easily extended to continuous spin 
models [131 ■ Nevertheless, the LL model and its vari- 
ants may be more closely connected to the Potts model 
than one may initially think. Note that for large p the 
Hamiltonian of Eq.(IT]) becomes increasingly Potts- like, 
in the sense that the pair interaction approaches a S- 
function: limp_>oo \di ■ dj\P = S{di,dj). This implies that 
neighbouring spins only interact when they are closely 
aligned and are otherwise indifferent to each other, just 
as in the Potts model. It has indeed been suggested that 
such models approximately resemble g-state Potts mod- 
els, with q ~ y/p |9|]. The observed trends in this work are 
certainly consistent with this interpretation. For all cases 
considered the strength of the transition increases with p, 
as manifested by the growing latent heats and interfacial 



tensions, exactly as in the Potts model with increasing q. 
Also the upward shift of e^o with p is consistent with the 
Potts model. However, it is clear that new theoretical 
approaches are needed to fully understand finite-size ef- 
fects at first-order transitions in the models studied here. 
We hope that the present simulation results may inspire 
such efforts. 
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